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SUMMARY 

Problem  4006  -  Cross  Correlation  Matrix  by  FFT 


The  problem  aim  was  to  produce  a  computer  program  for  determining 
cross  correlations  of  geosynchronous  satellite  images  by  means  of  a  two- 
dimensional  Fast  Fourier  Transform.  Features  to  be  addressed  in  deciding 
on  the  coding  for  the  program  were  accuracy  and  speed. 

A  computer  program  has  been  written  which  reads  data  from  a  perman¬ 
ent  file,  sets  up  arrays  of  specified  size,  and  at  the  user's  option  smooths 

2 

the  edges  of  the  tw-dimensional  array  via  a  cos  weighting  method.  The 
program  then  uses  our  two-dimensional  FFT  algorithm  to  compute  the  cross 
correlation  matrix  for  the  two  arrays.  An  option  is  available  which  searches 
by  row  and  column  to  locate  the  maximum  entry  in  the  row/ column  in  the  hope 
that  this  information  will  indicate  the  width  of  the  peaks. 

The  algorithm  used  in  our  routine  was  bench-marked  against  three 
commonly  used  FFT's.  In  the  most  competitive  case,  our  routine  was 
approximately  5  times  faster;  in  the  least  competitive,  we  were  approximate¬ 
ly  20  times  faster. 

Test  data  were  generated  and  run.  The  output  was  fully  consistent 
with  theoretical  analyses.  The  test  case  corresponded  to  a  two-dimensional 
square  wave.  This  was  chosen  for  its  simplicity  of  structure,  while  at 
the  same  time  providing  a  solid  check  on  the  accuracy  and  working  order  of 
the  algorithm. 

While  the  computer  program  was  being  developed,  alternative 
approaches  were  investigated.  Specific  attention  was  given  to  methods  which 
would  portentially  increase  processing  speed.  The  contemporary  journals  on 
information  and  image  processing  contain  an  abundance  of  references  to 
Walsh,  Haar  and  other  orthogonal  families  of  functions.  The  great  strength 
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of  the  Walsh  and  Haar  transforms  is  the  decrease  in  processing  times  by 
factors  of  2-20  over  comparable  times  using  the  FFT  algorithm.  After 
detailed  study  it  was  decided  that  these  methods  could  not  profitably  be 
applied  to  this  problem  at  this  time.  The  principal  drawback  is  the  state 
of  the  art  in  interpreting  and  decoding  the  convolution  and  deconvolution 
via  the  discrete  fast  Haar  and  fast  Walsh  transforms.  Furthermore,  since 
the  transform  variable  space  does  not  have  the  same  "intuitively  pleasing" 
relationship  as  do  time  and  frequency  (Fourier  transform  pairs) ,  it  is  not 
at  all  obvious  which  way  the  interpretations  should  (or  could)  be  pursued. 
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Let  f(x)  be  a  continuous  function  of  a  real  variable  x.  The 
Fourier  transform  of  f(x),  denoted  by  ?{f(x)|  ,  is  defined  by  the  equation 


3{f(x)  } 
where  j  * 


F(u)  f  - 
J- 


-J2nuxj 


f(x)  exp  - j  2nux  I  dx 


(1-1) 


Given  F(u) ,  f(x)  can  be  obtained  by  using  the  inverse  Fourier 
transform 


3  _1  |  F(u)  } 


f(x) 


r 


F(u)  exp  Jj2nux  du 


(1-2) 


Equations  (1-1)  and  (1-2) ,  which  are  called  the  Fourier  transform  pair,  can 
be  shown  to  exist  if  f(x)  is  continuous  and  integrable  and  F(u)  is  integrable. 
These  conditions  are  almost  always  satisfied  in  practice. 

We  will  be  concerned  with  functions  f(x)  which  are  real.  The 
Fourier  transform  of  a  real  function,  however,  is  generally  complex;  that  is. 


F(u) 


R(u)  +  jl(u) 


(1-3) 


where  R(u)  and  I(u)  are,  respectively,  the  real  and  imaginary  components  of 
F(u).  It  is  often  convenient  to  express  Eq.  (1-3)  inexponential  form: 


F(u)  -  | 

r(u)| 

(1-4) 

where 

j F(u)|  - 

^r200  +  i2(u)  j  1/2 

(1-5) 

and 

®(u) 

tm  [tw  ] 

(1-6) 
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The  magnitude  function  F(u)  is  called  the  Fourier  spectrum  of  f(x),  and  4  (u) 
its  phase  angle.  The  square  of  the  spectrum 

E(u)  -  j  F(u)  |  2  -  R2  (u)  +  I2  (u)  (1-7) 

is  commonly  referred  to  as  the  energy  spectrum  of  f(x). 

The  variable  u  appearing  in  the  Fourier  transform  is  often  called 
the  frequency  variable.  This  name  arises  from  the  fact  that,  using  Euler's 
formula,  the  exponential  term,  exp  [-j2rtux],  may  be  expressed  in  the  form: 

exp  [-2nux]  *  cos  2nux]  *  cos  2nux  -  j  sin  2nux  (1-8) 

>- 

If  we  interpret  the  integral  in  Eq.  (1-1)  as  a  limit-summation  of  discrete 
terms,  it  is  evident  that  F(u)  is  composed  of  an  infinite  sum  of  sine  and 
cosine  terms,  and  that  each  value  of  u  determines  the  frequency  of  its 
corresponding  sine-cosine  pair. 

For  the  step  function  shown  in  Fig  1  its  Fourier  transform  is  ob¬ 
tained  from  Eq.  (1-1)  as  follows: 

/■<x> 

f(x)  exp  [-j2rtux]  dx 

A  exp  I-j2nux]  dx 

sin  (rtuX)  e 
which  is  a  complex  function. 


F(u)  *  J 

-oo 

-/X 


A 

n  u 
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..... 


IH«)I 


(b) 


Figure  1.  A  step  function  and  its  Fourier  spectrum. 


The  Fourier  transform  can  he  easily  extended  to  a  function  f(x,y) 
of  two  variables.  If  f(x,y)  is  continuous  and  integrable,  and  F(u,v)  is 
integrable,  we  have  that  the  fillowing  Fourier  transform  pair  exists: 


and 


F  |f(x,y)|  -  F(u,v)  mJ^  f(*»y)  exp  I-j2n(ux  +  vy)]  dx  dy  (1-9) 

00 

F-1  F(u,v)  -  f(x,y)  -  J  J  F(u,v)  exp  [jtt(ux  +  vy)  ]  du  dv  (1-10) 


where  u  and  v  are  the  frequency  variables. 

As  in  the  one-dimensional  case,  the  Fourier  spectrum,  phase,  and 
energy  spectrum  are,  respectively,  given  by  the  relations: 


| F(u,v)  |  -  IR2  (u,v)  +  I2  (u,v)]  1/2 

(1-11) 

<b  (U,v)  -  tan  [r^] 

(1-12) 

and 

E(u,v)  -  R2  (u,v)  +  I2  (u,v) 

(1-13) 
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Example:  The  Fourier  transform  of  the  function  shown  in  Fig.  2(a)  is  given 

by: 


f(u,V)  -  jy 


—to 


AXY 


f(x,y)  exp  [-j2rc(ux  +  vy)]  dx  dy 

sin(nuX)  e  •^Tlu^  sin  (nvY)  e 

(nuX)  (tcvY) 


A*.y) 


A 


(a) 


l^v.  w)l 


I 


Figure  2(a). 

A  two-dimensional  function,  (b)  its  Fourier  spectrum. 


The  spectrum  is  given  by 


F(u,v) 


AXY 


sin  (ttuX) 
(nuX) 


sin  (TtvY) 
(tcvY) 


A  plot  of  this  is  shown  in  Fig.  2(b)  in  three-dimensional  perspective. 
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THE  DISCRETE  FOURIER  TRANSFORM 

Let  us  now  assume  that  a  continuous  function  f(x)  is  discretized  into 

a  sequence  Jf(x  ),  f(x  +  x) ,  f(x  +  2Ax) ,  ...  ,  f(x  +  [N  -  1]  Ax)  /  by 
C  o  o  o  o  j 

taking  N  samples  Ax  units  apart,  and  that  we  define 


f(x)  ■  f(x  +  xAx)  (2-1) 

o 

where  x  now  assumes  the  discrete  values  0,  1,  2,  ....  N  -  1.  In  other  words, 
the  sequency  |f(0),  f(l),  f(2),  ...  ,  f(N  -  1 )j  will  be  used  to  denote  any 
N  uniformly  spaced  samples  from  a  corresponding  continuous  function. 


With  the  above  notation  in  mind,  we  have  that  discrete  Fourier  trans¬ 
form  pair  that  applies  to  sampled  functions  is  given  by 


!  N"1 

F(u)  -  —  2  f(x)  exp  I  -j2Ttux/N] 

x=o 

for  u  =  0,  1,  2 . N  -  1,  and 

N-l 

f(x)  -  E  F(u)  exp  [  j2rcux/N] 
u=o 

for  x  =  0,  1,  2,  ...»  N-l. 


(2-2) 


(2-3) 


The  values  of  u  •  0,  1,  2,  ...,  N-l  in  the  discrete  Fourier  transform 
given  in  Eq.  (2-2)  correspond  to  samples  of  the  continuous  transform  at 
values  0,  Au,  2Au,  ....  (N-l)Au.  In  other  words  we  are  letting  F(u)  repre¬ 
sent  F(uAu) .  This  notation  is  similar  to  that  used  for  the  discrete  f(x), 
with  the  exception  that  the  samples  of  F(u)  start  at  the  origin  of  the  fre¬ 
quency  axis.  It  can  be  shown  that  Au  and  Au  are  related  by  the  expression 
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For  the  two  variable  case  the  discrete  Fourier  transform  pair  is 
given  by  the  equations 


F<u>v)  • 


M-l 

£ 

x-o 


N-l 

£  f(x  y)  exp  [~j2n(ux/M  +  vy/N)) 
yo 


(2-5) 


for  u  »  0,  1,  2,  ...M 

M-l  N-l 
f(x,y)  =  £  £ 

u=o  v«o 


1 ,  v  ■  0 ,  1,  2 ,  ...  ,  N  -  1 ,  and 

F(u,v)  exp  [j2n(ux/M  +  vy?N)) 


for  x  •  0,  1,  2, 


N-l  and  y  -  0 ,  1 ,  2 ,  ...N-l. 


(2-6) 


Sampling  of  a  continuous  function  is  now  in  a  two-dimensional  grid 
with  divisions  of  width  Ax  and  Ay  in  the  x  and  y  axis,  respectively.  As  in 
the  on-dimensional  case,  the  discrete  function  f(x,y)  represents  samples  of 
the  function  f(xQ  +  xAx,  yQ  +  yAy)  for  x**0,  1,  2,  ...  ,  M  -  1  and  y  - 
0,  1,  2,  ...  ,  N  -  1.  Similar  comments  hold  for  F(u,v).  The  sampling 
increments  in  the  spatial  and  frequency  domains  are  related  by 


and 


1 

MAx 


1 

NAy 


(2-7) 

(2-8) 


When  images  are  sampled  in  a  square  array  M  ■  N  and 
N-l  N-l 

F(u,v)  =  ^  £  E  f(xy)  exp  [-j2n(ux  +  vy)  /  N] 

x=o  y=*o 


(2.9) 


for  u,  v  =  0,  1,  2,  ...  ,N-1  and 
n  NhI  N-l 

£  £  F(u,v)  exp  [j2n(ux  +  vy)  /  N] 


f(x,y)  =  n 


(2.10) 


u=o  v=o 


for  x,  y  =  0,  1,  2,  ...,  N  -1.  Note  that  in  the  case  we  have  included  a  1/N 
term  in  both  expressions.  Since  F(u,v)  and  f(x,y)  are  a  Fourier  transform 
pair,  the  grouping  of  these  constant  multiplicative  terms  is  arbitrary. 
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EXAMPLE:  As  an  illustration  of  the  use  of  Eqs.  (2-2)  and  (2-3),  consider 
the  functions  shown  in  Fig.  5(a).  If  this  function  is  sampled  at  the  argu¬ 
ment  values  x  *  0.5,  x,  ■  0.75,  x,  *  1.0,  x,  ■  1.25,  and  if  the  argument  is 

O  1  4  J 

redefined  as  discussed  above,  we  obtain  the  discrete  function  shown  in 
Fig.  3(b). 


Figure  3. 

A  simple  function  and  samples  in  the  x  domain. 
In  (a)  x  is  a  continuous  variable;  in  (b)  x  is 
discrete. 


Application  of  Eq.  (2-2)  yields  the  following: 

1  3 

F(0)  -  E  f(x)  exp  [0] 

x*o 

-  ^  f(0)  +  f(l)  +  f (2)  +  f (3)  -  3.25 

1  3 

F(l)  -  j  E  f  (x)  exp  [-j  2txx/4] 

x*o 

-  y  [2e°  +  3e  ~^2  +  4e  "jTI  +  4e  “J3n/2^  -  i  [-2+j]. 

4  ** 


Notice  that  all  values  of  f(x)  contribute  to  each  of  the  four  terms 
of  the  discrete  Fourier  transform.  Conversely,  all  terms  of  the  transform 
contribute  in  forming  the  inverse  transform. 


SOME  PROPERTIES  OF  THE  TWO-DIMENSIONAL  FOURIER  TRANSFORM 

While  our  primary  interest  is  in  two-dimensional,  discrete  transforms, 
the  underlying  concepts  of  some  of  these  properties  are  much  easier  to  grasp 
is  they  are  presented  in  their  one-dimensional,  continuous  form. 

PERIDIOCITY  AND  CONJUGATE  SYMMETRY 


The  discrete  Fourier  transform  and  its  inverse  are  periodic  with 
period  N;  that  is, 

F(u,v)  =  F(u  +  N,  v)  “  F(u,v  +  N)  =  F(u  +  N,  v  +  N)  (3.1) 

Although  Eq.  (3-1)  points  out  that  F(u,v)  repeats  itself  for  an  infinite 
number  of  values  of  u  and  v,  only  the  N  values  of  each  variable  in  any 
one  period  are  required  to  obtain  f(x,  y)  from  F  (u,  v) .  In  other  words, 
only  one  period  of  the  transform  is  necessary  to  completely  specify  F(u,v) 
in  the  frequency  domain.  Similar  comments  hold  for  f(x,y)  in  the  spatial 
domain . 


The  Fourier  transform  also  exhibits  conjugate  symmetry  since 


F  (u,v)  =  F*  (  -u,  -v) 


(3-2) 


or  more  interestingly, 

|F(u,v)  |  =  F  |  (-u,  -v)  | 


(3-3) 


As  mentioned  earlier,  it  is  often  of  interest  to  display  the  magnitude  of  the 
Fourier  transform  for  interpretation  purposes. 
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The  same  observations  hold  for  the  magnitude  of  the  two-dimensional  i 

Fourier  transforms,  with  the  exception  that  the  results  are  considerably  ' 

more  difficult  to  interpret  if  the  origin  of  the  transform  is  not  shifted  to 
the  frequency  point  (N/2,  N/2). 

D1STRIBUTIVITY  AND  SCALING 


It  follows  directly  from  the  definition  of  the  continuous  or  discrete 
transform  pair  that 

F  ^(x.y)  +  f2  (x,y)J  “  F  ^  (x,y)J  +  F  |f2  (x,y)j  and,  in  general 

that 

F  (x.y)  .  f2  (x,y)  J4  F^f][  (x,y)J  .  F  £f2  (x,y)j 

In  other  words,  the  Fourier  transform  and  its  inverse  are  distributive  over 
addition,  but  not  over  multiplication. 

It  Is  also  easy  to  show  that  for  two  scalars  a  and  b, 

af (x.y)  aF(u,v) 

and 

f(ax,by)  F  *u/a,  v/b) 

The  two-dimensional,  discrete  convolution  is  formulated  by  letting 
f(x,y)  and  g(x,y)  be  discrete  arrays  of  size  A  X  B  and  C  X  D,  respectively. 
As  in  the  one-dimensional  case,  these  arrays  must  be  assumed  periodic  with 
some  period  M  and  N  in  the  x  and  y  directions,  respectively.  Wraparound 
error  in  the  individual  convolution  periods  is  avoided  by  choosing 

M>/A  +  C  -1  (3-9) 

and 

N>/B  +  D  -  1  (3.10) 
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The  periodic  sequences  are  formed  by  extending  f(x,y)  and  g(x,y)  as 
follows: 


fe(x.y) 


f  (x,y)  0  £x  £A  -1  and  0,$  y  £  B  -  1 


As<  x£M  -  1  or 


B$y$N  -  1 


and 


ge(x,y) 


f(x,y)  0^x<C  -  1  and  0^y{D  -  1 


CN<x/ M  -  1  or  D{y(N  -  1 


The  two-dimensional  convolution  of  f£(x,y)  and  g£(x,y)  is  given  by 
the  relation 


fe(x,y)*ge(x,y) 


M-l 

E 

m-0 


N-l 

E  f  (m,n)  g  (x  -  m,  y  -  n)  (3.11) 
n-0  e  e 


for  x  =  0,  1,  2,  ...  ,M-1  and  y  ■  0,  1,  2,  ...N-l.  The  M  X  N  array 
given  by  this  equation  is  one  period  of  the  discrete,  two-dimensional  convo¬ 
lution.  If  M  and  N  are  chosen  as  Indicated  above  this  array  is  guaranteed 
to  be  free  of  interference  from  other  adjacent  periods.  All  computations 
use  these  extended  functions  fe(x,y)  and  gg(x,y). 

AVERAGE  VALUE 


A  widely-used  definition  of  the  average  value  of  a  two-dimensional 
discrete  function  is  given  by  he  expression 

N-l  N-l 

f(x,y)  m  E  E  f(x,y)  -  |  f(0,0)  (3.4) 

N  x-0  y-0 
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CONVOLUTION  and  CORRELATION 


We  will  now  outline  two  Fourier  transform  relationships  which  consti¬ 
tute  a  fundamental  link  between  the  spatial  and  frequency  domains.  These 
relationships,  called  convolution  and  correlation,  are  of  fundamental 
importance  in  developing  a  firm  understanding  of  image  processing  techniques 
based  on  the  Fourier  transform. 

CONVOLUTION 


The  convolution  of  two  functions  f(x)  and  g  (x) ,  denoted  by  f(x)*g(x), 
is  defined  by  the  integral 


■P 

^  — on 


f(x)*g(x)  -  I  f  (  a)  g(x  -a)  da 
-oo 

where  a  is  the  variable  of  integration. 


(3-5) 


Two-dimensional  convolution  is  analogous  in  form  to  Eq.  (3-5).  For 
two  functions  f(x,y)  and  g(x,y) ,  we  have 


f(x,y)*g(x,y) 


■//. 


f(a,3  )g(x-a,y-fi)  dodp  (3-6) 


The  convolution  theorem  in  two-dimensions  is  summarized  by  the  relations 


f(x»y)*g(x,y)  F(u,v)G(u,v) 

and 

f(x.y)g(x,y)  F(u,v) *G(u,v) 


(3-7) 


(3-8) 


From  a  practical  point  of  view,  it  is  often  more  efficient  to  compute 
the  discrete  convolution  in  the  frequency  domain  instead  of  using  Eq.  (3-11) 
directly.  The  procedure  is  to  compute  the  Fourier  transforms  of  fg(x,y)  and 
gg(x,y)  byusing  a  fast  Fourier  transform  (FFT)  algorithm.  The  two  transforms 
are  then  multiplied  and  the  inverse  Fourier  transform  of  the  product  will 
yield  the  convolution  function.  A  comparison  by  Brigham  (1974)  shows  that, 
for  one-dimensional  arrays,  the  FFT  approach  is  faster  if  the  number  of  points 
is  greater  than  32.  Although  this  figure  is  dependent  on  the  particular 
machine  and  algorithms  used,  it  is  well  below  the  number  of  points  in  a  row 
or  column  of  a  typical  image. 

CORRELATION 


The  correlation  of  two  continuous  functions  f(x)  and  g(x)  denoted  by 
f(x)  og  (x) ,  is  defined  by  the  relation 

f(x)°g(x)  ■  f  (  o.)  g  (x  +  a)  da  (3-12) 


To  perform  correlation  one  slides  g(x)  by  f(x)  and  integrate  the  product 
from  -oo  to  oo  for  each  displacement  x. 

The  discrete  equivalent  of  continuous  correlation  is  defined  by 
M-l 

f  (x)°ge(x)  «  S  fg  (m)  gg  (x  +  m)  (3-13) 

'  m“0 


for  x  *  0,  1,  2,  ...,  M  -  1.  The  same  comments  made  above  with 

regard  to  f  (x)  and  g  (x) ,  the  assumed  periodicity  of  these  functions,  and 
e  e 

the  choice  of  values  for  M,  apply  to  this  case. 


Similar  expressions  hold  for  two  dimensions.  If  f(x,y)  and  g(x,y) 
are  functions  of  continuous  variables,  their  correlation  is  defined  as 


*(xjy)°  g  (x,y) 


Jf 

-oo 


f(a£)  g  (x  +  a,  y  +£  ) 


d  a  d(3 


(3-14) 


v 
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and  for  the  discrete  case  we  have 


M-l  N-l 

f  (x,y)°g  (x,y)  -  2  E  f  (m,  n)  g  (x+m,  y+n)  (3-15) 

e  m-0  n-0  e  e 

for  x  -  0,  1,  2,  ...  ,  M  -  1  and  y  ■  0,  1,  2,  ...»  N  -  1.  As  in  the  case  of 
discrete  convolution,  fg  (x,y)  and  gg  (x,y)  are  extended  functions,  and  M  and 
N  are  chosen  according  to  avoid  wraparound  error  in  tie  periods  of  the  correla¬ 
tion  function. 

For  both  the  continuous  and  discrete  cases  the  following  correlation 
theorem  holds: 

f  (x,y)°  g  (x,y)  ^  F(u,v)  G*  (u,v) 

and 

f(x,y)  g*  (x,y)  F(u,v)°  G  (u,v) 

where  denotes  the  complex  conjugate.  Note:  for  discrete  variables,  all 
functions  are  assumed  to  be  extended  and  periodic. 

One  of  the  principal  applications  of  correlation  in  image  processing 
is  in  the  area  of  template  or  prototype  matching,  where  the  problem  is  to 
bind  the  closest  match  between  a  given  unknown  image  and  a  set  of  images  of 
known  origin.  One  approach  to  this  problem  is  to  compute  the  correlation 
between  the  unknown  and  each  of  the  known  images.  The  closest  match  can 
then  be  found  be  selecting  the  image  that  yields  the  correlation  function 
with  the  largest  value.  Since  the  resultant  correlations  are  two-dimensional 
functions,  this  involves  searching  for  the  largest  amplitude  of  each  function. 
As  in  the  case  of  discrete  convolution,  the  computation  of  f  (x,y)°ge  (x,y) 
is  often  more  efficiently  carried  out  in  the  frequency  domain  using  an  FFT 
algorithm  to  obtain  the  forward  and  inverse  transforms. 

THE  FAST  FOURIER  TRANSFORM 

The  number  of  complex  multiplications  and  additions  required  to  imple- 

2 

ment  Eq.  (2-2)  is  proportional  to  N  .  This  can  be  seen  easily  by  noting 
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that,  for  each  of  the  N  values  of  u,  expansion  of  the  summation  sign  re¬ 
quired  N  complex  multiplications  of  f (x)  hy  exp  [-.)2n.ux/N]  can  be  computed 
once  and  siored  in  a  table  for  all  subsequent  applications.  For  this 
reason,  the  multiplication  of  u  by  x  in  these  terms  is  usually  not  considered 
a  direct  part  of  the  implementation. 

In  this  section  it  is  shown  that,  by  properly  decomposing  the  number 
of  multiply  and  add  operations  it  can  be  made  proportional  to  N  log2  N.  The 

decomposition  procedure  is  called  the  Fast  Fourier  transform  (FFT)  algorithm. 

2 

The  reduction  in  proportionality  from  N  to  N  log^  N  multiply/add  operations 
represents  a  significant  savings  in  computation  effort,  as  shown  by  the 
figures  in  the  table.  It  is  evident  from  this  table  that  the  FFT  approach 
offers  a  considerable  computational  advantage  over  a  direct  implementation 
of  the  Fourier  transform,  particularly  when  N  is  relatively  large.  For 
example,  a  direct  implementation  of  the  transform  for  N-8192  requires  on  the 
order  of  three-quarters  of  an  hour  in  a  machine  such  as  an  IBM  7094.  By 
contrast,  the  same  job  can  be  done  in  this  machine  in  about  5  sec.  using  an 
FFT  algorithm. 


N 

N2 

(Conventional  FT) 

N  log  N 
(FFT) 

Computational  Advantage 
(N/log2N) 

2 

4 

2 

2.00 

4 

16 

8 

2.00 

8 

64 

24 

2.67 

16 

256 

64 

4.00 

32 

1,024 

160 

6.40 

64 

4,096 

384 

10.67 

128 

16,384 

896 

18.29 

256 

65,536 

2,048 

32.00 

512 

262,144 

4,608 

56.89 

1024 

1,048,576 

10,240 

102.40 

2048 

4,194,304 

22,528 

186.18 

4096 

16,777,216 

49,152 

341.33 

8192 

67,108,864 

106,496 

630.15 

As  pointed  out  previously,  a  two-dimensional  Fourier  transform  can  be 
readily  computed  by  a  series  of  applications  of  the  one-dimensional  transform. 


THE  INVERSE  FFT 

Thus  far,  little  has  been  said  concerning  the  inverse  Fourier  transform. 
It  turns  out  that  any  algorithm  for  implementing  the  discrete  forward 
transform  can  also  be  used  (with  minor  modifications  in  the  input)  to  compute 
the  inverse. 

For  a  two-dimensional  square  array  we  take  the  complex  conjugate,  that 
is. 


f*(x,y)  -  ^  Nj1  Nj1  F*(u,v)  exp  [-j2n(ux  +  vy)/N] 

u-o  v“0 

which  we  see  is  in  the  form  of  the  two-dimensional  forward  transform.  It 
follows,  therefore,  that  if  we  input  F*  (u,v)  into  an  algorithm  designed  to 
compute  the  forward  transform,  the  result  will  be  f*  (x,y).  By  taking  the 
complex  conjugate  of  this  result  we  obtain  f(x,y).  In  the  case  where  f(x) 
or  f(x,y)  are  real,  the  complex  conjugate  operation  is  unnecessary  since 
f(x)  *  f*(x)  and  f(x,y)  ■  f*(x,y)  for  real  functions. 

The  fact  that  the  two-dimensional  transform  is  usually  computed  by 
successive  passes  of  the  one-dimensional  transform  is  a  frequent  source  of 
confusion  when  using  the  above  technique  to  obtain  the  inverse.  When  using 
a  one-dimensional  algorithm  to  compute  the  two-dimensional  inverse,  the  method 
is  not  to  compute  the  complex  conjugate  after  each  row  or  column  is  processed. 
Instead,  the  function  F*(u,v)  is  treated  as  if  it  were  f(x,y)  in  the  forward, 
two-dimensional  transform  procedure.  The  complex  conjugate  of  the  result 
(if  necessary)  will  yield  the  proper  inverse  f(x,y). 
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A  number  of  methods  exist  which  generalize  the  Fourier  transform 
either  by  considering  the  class  of  Fast  Unitary  Transforms  (which  include 


the  FFT)  or  by  considering  group  characters.  The  practical  interest  here 
is  in  the  computational  efficiency  inherent  in  these  more  general  transforms. 

For  an  arbitrary  sequence  of  functions,  the  Gram-Schmidt  process 
generates  a  sequence  of  orthogonal  functions.  Any  continuous  function  can 
be  expressed  via  this  orthogonal  set  with  minimum  error;  for  certain 
classes  of  functions  the  convergence  is  uniform.  Most  inportant  among  these 
are  the  Baar  functions,  which  assume  2  values,  and  the  Walsh  functions  with 
values  +  1.  Both  Both  have  discrete  analogies  and  Fast  Discrete  Haar/Walsh 
transforms  are  computable,  FHT  and  FOT.  Because  of  the  simplicity  of  the" 
basic  functions  much  greater  computational  savings  can  be  obtained  than  from 
FFT  (up  to  30  times  faster  than  FFT)  . 

Advantageous  use  can  be  made  of  FHT/FWT  in  certain  applications;  e.g. 
in  data  transmission/reconstruction  in  which  one  represents  a  given  signal 
in  some  sense  and  reconstruct  it  from  the  minimal  representation.  A  number 
of  serious  difficulties  arise  with  these  transforms  due  to  the  fact  that 
the  relation  with  the  circle  has  been  lost.  For  example:  1)  no  natural 
interpretation  in  terms  of  frequency  exists  (the  "sequency"  viewpoint  of 
Harmuth  for  Walsh  functions  lacks  physical  meaning);  2)  due  to  absence  of 
the  circle  relationship  the  important  convolution  theorem  is  not  available 
(forced  analogies  to  a  convolution  theorem  via  dyadic  convolutions  have  been 
made  but  their  interpretations  are  not  clear) . 

The  striking  advantages  of  FWT  and  FHT  over  the  usual  FFT  in  computa¬ 
tional  effort  should  motivate  further  investigation  in  this  area. 

The  historical  situation  regarding  orthogonal  functions  at  the 
beginning  of  the  twentieth  century  was  one  of  well  known  and  useful  kinds  of 
such  functions:  the  trigonometric  functions  which  occur  in  Fourier  series; 
orthogonal  polynomials  such  as  those  of  Legendre,  Hermite,  and  Laguerre; 


Bessel”s  functions,  the  Sturm-Liouville  series,  and  other  special  functions. 
But,  there  was  no  general  theory  embracing  all  such  systems  of  functions. 

The  Hungarian  mathematician,  Alfred  Haar,  was  concerned  with  con¬ 
vergence  properties  of  series  of  orthogonal  functions,  and  also  constructing 
a  new  set  (now  called  the  Haar  system)  of  such  functions.  He  defined  a  set 
of  orthogonal  functions  each  taking  essentially  only  two  values  such  that 
the  formal  expansion  of  an  arbitrary  continuous  function  in  those  functions 
converges  uniformly  to  the  given  function,  a  property  not  possessed  by 
orthogonal  sets  known  up  to  that  time. 

..  In  1923,  J.L.  Walsh  published  a  set  of  orthogonal  functions  which 

are  complete  on  the  interval  [0.1];  they  take  only  the  values  +  1,  and  are 
similar  in  oscillation  and  many  other  properties  to  the  trigonometric 
functions.  They  have  turned  out  to  have  important  practical  applications 
in  calculation. 

The  limits  of  the  usefulness  of  these  functions  both  in  theoretical 
work  and  in  engineering  applications  still  seem  to  be  undetermined. 

Traditionally,  the  theory  of  communication  has  been  based  on  the 
complete,  orthogonal  system  of  sine  and  cosine  functions.  The  concept  of 
frequency  is  defined  as  the  parameter  f  in  sin  2n  ft.  and  cos  2ti  ft.  The 
question  arises  whether  there  are  other  systems  of  functions  on  which 
theories  of  similar  scope  can  be  based,  and  that  lead  to  equipment  of 
practical  interest. 

The  parameter  in  ^2sin  2in9  and  /2~cos  2in9  gives  the  number  of 
oscillations  in  the  interval  -1/2  <9  <1/2  (that  is,  the  normalized  frequency 
i=fT).  One  may  interpret  i  as  "one  half  the  number  of  zero  crossings  per 
unit  time"  rather  than  as  "oscillations  per  unit  time".  (The  zero  crossing 
at  the  left  side,  9  =  -1/2,  but  not  the  one  at  the  right  side,  9  =  +1/2,  of 
the  time  interval  is  counted  for  sine  functions). 
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The  parameter  i  also  equals  one  half  the  number  of  zero  crossings  in 
the  interval  -1/2  £  9  ^ 1/2  for  Walsh  functions.  In  contrast  to  sine-cosine 
functions,  the  sign  changes  are  not  equidistant.  If  i  is  not  an  integer, 
then  it  equals  "one  half  the  average  number  of  zero  crossings  per  unit  time". 
The  term  "normalized  sequency"  has  been  introduced  for  L  and  ©  ■  i/T  is 
called  the  normalized  sequency.  Sequency  in  zps  »  1/2  (average  number  of 
zero  crossings  per  second) . 

The  general  form  of  a  sine  function  V  sin  (2nft+c£>  contains  the 
parameters  amplitude  V,  frequency  f,  and  phase  angle  a  The  general  form 
of  a  Walsh  function  V  sal  (®T,  t/T  +  to/T)  contains  the  parameters  amplitude, 
V,  sequency,  ®  ,  the  delay,  tQ,  and'  time  base,  T.  The  normalized  delay, 
t^/T,  corresponds  to  the  phase  angle.  The  time  base,  T,  is  an  additional 
parameter  and  it  causes  a  major  part  of  the  differences  in  the  applications 
of  sine-cosine  and  Walsh  functions. 


So  far,  Walsh  functions  are  the  only  known  functions  with  desirable 
features  comparable  to  sine-cosine  functions  for  use  in  communications. 
Development  of  semi-conductor  technology  has  imparted  practical  interest 
in  them  at  this  time.  Generally* speaking,  the  transition  from  sine-cosine 
functions  to  other  complete  systems  means  a  transition  from  linear,  time- 
invariant  components  and  equipment  to  linear,  time-variable  components  and 
equipment,  which,  of  course,  constitute  a  much  larger  class.  The  mathemati¬ 
cal  theory  of  Walsh-Fourier  analysis  corresponds  to  the  Fourier  analysis 
used  for  sine-cosine  functions.  There  is  no  theory  of  similar  scope  for 
block  pulses,  because  they  are  incomplete. 


The  sal  and  cal  transforms  of  Walsh-Fourier  analysis  are  defined  by 


a  (U) 


a  (u) 


s: 

/: 


F  (9)  sal  (u,0)  d© 


F  (0)  cal  (u,©)  d© 
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a)  Orthonormal  systems  of  functions. 
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b)  List  of  features  and  applications  of  sine-cosine  functions,  Walsh 
functions,  and  block 
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Walsh-Function  Filters  -  For  a  sequency  low  pass  filter  based  on  Walsh 
functions  the  input  signal,  F(9),  is  transformed  into  a  step  function. 

Ff T(9) ,  with  steps  of  a  certain  width,  by  integrating  F(9)  during  an  interval 
equal  to  the  step  width.  The  amplitudes  of  the  steps  are  chosen  so  that 
Fft(9)  yields  a  least-mean-square  approximation  of  F(9) .  In  addition, 

Ftt(9)  is  delayed  with  respect  to  F(9)  by  one  step  width. 

The  number  of  samples  obtained  ia  equal  to  twice  the  cut-off 
sequency.  Hence,  the  sampling  theorems  of  Fourier  analysis  permit  the 
comparison  of  frequency  and  sequency  filters. 

-Theorems  for  the  multiplications  of  Walsh  functions  have  been 
proven.  These  are: 

cal  (k,  9)  cal  (1,  0)  «  cal |  k  i,  9l 

sal  (k,  9)  cal  (i,  9)  -  sal  |  i  &  (k-1)  +  1,9| 

sal  (k,  9)  sal  (i,  0)  -  cal  |(k-l)  Q  (i-1) ,  0| 

where  the  symbol  ©  indicates  modulo  2  addition.  Note  that  the  product  of 
two  Walsh  functions  yields  only  one  Walsh  function.  Therefore,  the  amplitude 
modulation  of  a  Walsh  carrier  yields  only  one  sequency  sideband  as  compared 
to  the  two  sidebands  obtained  when  a  sine  carrier  is  modulated.  A  typical 
application  of  the  multiplication  theorems  of  Walsh  functions  is  in  the 
design  of  sequency-bandpass  filters. 

Digital  Filtering  and  Multiplexing  -  One  of  the  most  promising  aspects  of 
Walsh  functions  is  the  case  with  which  filters  and  multiplex  equipment  can 
be  implemented  as  digital  circuits.  The  reason  is  that  numerical  Walsh- 
Fourier  transformation  and  numerical  sequency  shifting  of  signals  require 
summations  and  subtractions  only.  In  the  case  of  sine-cosine  functions,  the 
corresponding  operations  require  multiplications  with  irrational  numbers. 


A  digital  filter  based  on  Walsh  functions  can  be  readily  obtained. 
The  input  signal  passes  first  through  a  sequency  low-pass  filter  then 
transforms  it  into  a  step  function.  This  step  function  is  sampled  and  the 
samples  are  transformed  into  numbers  by  an  analog/digital  converter.  A 
series  of  these  numbers  is  stored  in  a  digital  storage.  A  Walsh-Fourier 
transform  of  this  series  is  obtained  by  performing  certain  additions  and 
subtractions  in  an  arithmetic  unit.  Some  or  all  of  the  obtained  co¬ 
efficients,  that  represent  sequency  components,  may  be  suppressed  or 
altered  -  in  effect,  a  filtering  process.  An  inverse  Walsh-Fourier  trans¬ 
form  yields  the  filtered  signal  as  a  series  signal  by  digital/analog  con¬ 
verter.  Since  there  is  a  fast  Walsh-Fourier  transform  just  as  ther  is  a 

-  >■ 
fas  Fourier  transform,  the  arithmetic  operations  in  a  digital  sequency  filter 

are  not  only  simpler  than  in  a  digital  frequency  but  can  be  performed  faster. 

One  of  the  features  of  Walsh  functions  that  makes  them  of  some 
interest  in  signal  processing  is  the  fact  that  their  amplitudes  are  given 
precisely  by  a  single  bit,  so  that  their  use  does  not  directly  contribute  to 
roundoff  noise.  The  basis  vectors  of  symmetry  analysis  offer  the  same 
attraction  with,  additionally,  for  low  orders  of  input  data  frames  N,  some 
economy  of  computations  by  reason  of  the  zeros. 

OTHER  IMAGE  TRANSFORMS 


The  Fourier  transform  is  the  transform  most  often  used  in  image 
processing  applications;  there  are  other  transforms  which  are  also  of 
interest  in  this  area. 

The  one-dimensional,  discrete  Fourier  transform  is  one  of  a  class 
of  important  transforms  which  can  be  expressed  in  terms  of  the  general 
relation 

N-l 

T(u)  =  E  f<x)g(*,u)  (5-1) 

x*0 
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where  T(u)  is  the  transform  of  f(x),  g(x,u)  ie  the  forward  transformation 
kernel,  and  u  assumes  values  in  the  range  0,  1,  N-l.  Similarly,  the 
inverse  transform  is  given  by  the  relation 

N-l 

f(x)  -  £  T(u)h(x,u)  (5-2) 
u-0 

where  h(x,u)  is  the  inverse  transformation  kernel  and  x  assumes  values  in 
the  ranges  0,  1,  ....  N-l.  The  nature  of  a  transform  is  determined  by  the 
properties  of  its  transformation  kernel. 


For  two  dimensional  square  arrays  the  forward  and  inverse  transforms 
are  given  by  the  equations 


and 


T(u,v) 


N-l  N-l 

£  £  f (x,y)g(x,y,u,v) 

x-0  y-0 


(5-3) 


f(x,y) 


N-l  N-l 

£  £  T(u,v,)h(x,y,u,v) 

u-0  v**0 


(5-4) 


where,  as  above,  g(x,y,u,v)  and  h(x,y,u,v)  are  called  the  forward  and 
inverse  transformation  kernels,  respectively. 


The  two  dimensional  Fourier  transform  has  the  kernel 


g(x,y,u,v)  -  ^  exp  [~j2n  (ux  +  vy)  /nJ 


which  is  separable  and  symmetric  since 


g(x,y,u,v)  -  g1(x,u)g1(y,v) 


j-mr  exp  iT-j2nux/Nj  —  ^=.  exp  [*-j2Ttvy/Nl 
/  N  /  N  L 
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It  is  easily  shown  that  the  inverse  Fourier  kernel  is  also  separable  and 
symmetric. 

A  transform  with  a  separable  kernel  can  be  computed  in  two  steps, 
each  requiring  a  one  dimensional  transform.  First,  the  one  dimensional 
transform  is  taken  along  each  row  of  f(x,y),  yielding 

N-l 

T(x,v)  -  S  f(x,y)g,(y,v)  (5-5) 

y-o 

for  x,v  “0,  1,  2,  ....  N-l.  Next,  the  one  dimensional  transform  is  taken 

along  each  column  of  T(x,v) ;  this  results  in  the  expression 

► 

N-l 

T(u,v)  *  2  T(x,v)g  (x,u)  (5-6) 

x*0 

for  u , v  =  0,  1,  2,  ....  N-l.  The  same  final  results  are  obtained  if  the 
transform  is  taken  first  along  each  column  of  f(x,y)  to  obtain  T(y,u)  and 
then  along  each  row  of  the  latter  function  to  obtain  T(u,v).  Similar 
comments  hold  for  the  inverse  transform  if  h(x,y,u,v)  is  separable. 

If  the  kernel  g(x,y,u,v)  is  separable  and  symmetric,  Eq.  (5-3)  can 
also  be  expressed  in  the  following  matrix  form: 

T  =  AFA  (5-7) 

where  F  is  the  NXN  image  matrix,  A  is  an  NXN  symmetric  transformation  matrix 
with  elements  a^  =  g^  (i,j),  and  T  is  the  resulting  NXN  transform  for 
values  of  u  and  v  in  the  range,  0,  1,  2,  ...»  N-l. 

To  obtain  the  inverse  transform  we  pre-multiply  and  post-multiply 
Eq.  (5-7)  by  an  inverse  transformation  matrix  B. 


If  B  *  A  \  it  then  follows  that 
F  *  BTB 

which  indicates  that  the  digital  image  F  can  be  recovered  completely  from 
its  transform.  If  B  is  not  equal  to  A  \  then  we  obtain  an  approximation 
to  F,  given  by  the  relation 

F  -  BAFAB 

A  number  of  transforms,  including  the  Fourier,  Walsh,  and  Haar 
transforms,  can  be  expressed  in  this  form.  An  important  property  of  the 
resulting  transformation  matrices  is  that  they  can  be  decomposed  into 
products  of  matrices  with  fewer  non-zero  entries  than  the  original  matrix. 
This  result,  first  formulated  by  Good  (1958)  for  the  Fourier  transform, 
reduces  redundancy  and,  consequently,  the  number  of  operations  required  to 
implement  a  two-dimensional  transform.  The  degree  of  reduction  is  equiva¬ 
lent  to  that  achieved  by  an  FFT  algorithm,  being  on  the  order  of  N  log2  N 
multiply/add  operations  for  each  row  or  column  of  an  NXN  image. 

Walsh  Transform 


When  N  *  2n,  the  discrete  Walsh  transform  of  a  function  f(x), 

denoted  by  W(u),  is  obtained  by  substituting  the  kernel 


g(x,u)  =  -  II  (-1)  b1(x)bn_i_i^u^  (5-8) 


into  Eq.  (5-1).  In  other  words, 

,  N-l  N-l 

W(u)  -  £  E  f(x)  (-1)  b  (x)b  .(u)  (5-9) 

N  x-0  i-0  1  n"1"1 
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where  bfc(z)  is  the  kth  bit  in  the  binary  representation  of  z.  For  example, 
if  n  -  3  and  z  -  6  (110  in  binary),  we  have  that  bQ(z)  =  0,b^(z)  *  1,  and 
b2(z)  -  1. 

The  values  of  g(x,u),  excluding  the  l/N  constant  term,  are  listed 
below  for  N  *  8.  The  array  formed  by  the  Walsh  transformation  kernel  is 

Values  of  the  Walsh  trans¬ 
forms /tionkernel forN“8 . 


0  1  2  3  4  5  6  7 


0 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

1 

+ 

+ 

+ 

+ 

- 

- 

- 

- 

2 

+ 

+ 

- 

- 

+ 

+ 

- 

- 

3 

+ 

+ 

- 

- 

- 

- 

+ 

+ 

4 

+ 

- 

+ 

- 

♦ 

— 

+ 

- 

S 

+ 

- 

+ 

- 

- 

+ 

- 

+ 

6 

+ 

- 

- 

+ 

♦ 

- 

- 

* 

7 

+ 

- 

- 

+ 

- 

■f 

♦ 

- 

a  symmetric  matrix  whose  rows  and  columns  are  orthogonal.  These  properties, 
which  hold  in  general,  lead  to  an  inverse  kernel  which  is  identical  to  the 
forward  kernel,  except  for  a  constant  multiplicative  factor  of  l/N.  Thus, 
the  inverse  Walsh  transform  is  given  by 

N-l  n-1 

f(x)  *»  £  W(u)  n  (-1)  b  (x)b  (u)  (5-10) 

u=0  1=0 

Notice  that,  unlike  the  Fourier  transform  which  is  based  on  trigono¬ 
metric  terms,  the  Walsh  transform  consists  of  a  series  expansion  of  basis 
functions  whose  values  are  either  plus  or  minus  one. 

It  is  also  of  interest  to  note  that  the  forward  and  inverse  Walsh 
transforms  differ  only  by  the  l/N  term.  Thus,  any  algorithm  for  computing 
the  forward  transform  can  be  used  directly  to  obtain  the  inverse  transform 
simply  by  multiplying  the  result  of  the  algorithm  by  N. 
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The  forward  and  inverse  Walsh  transforms  are  also  equal  given  by 


W(u,v) 


and 


f (x,y) 


,  N-l  N-l 
J  2  E  f(x,y) 

x-0  y-0 


n-l 

n  (-1)  VX>bn-l-i(u)  +  bi(y)bn-l-i(v) 

i-0 

(5-11) 


N-l  N-l 

i  S  E  W(u,v) 
u-0  v=0 


n-l 

n 

i-0 


(-1)  v1(*)bn^1(u)rt1(y)bn.1_i(v) 

(5-12) 


Thus,  any  algorithm  which  is  used  to  compute  the  two  dimensional  forward  ' 
Walsh  transform  can  also  be  used  without  modification  to  compute  the  inverse 
transform. 


The  Walsh  transform  can  be  computed  by  a  fast  algorithm  identical  in 
form  to  the  successive  doubling  method  for  the  FFT.  The  only  difference  is 
that  all  exponential  terms  WN  are  set  equal  to  one  in  the  case  of  the  fast 
Walsh  transform  (FWT) . 

The  Walsh  transform  is  real,  thus  requiring  less  computer  storage 
for  a  given  problem  than  the  Fourier  transform,  which  is  in  general  complex 
valued. 
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INTRODUCTION 


The  program  which  controls  all  of  the  calculations  and  output  related 
to  this  task  is  called  WETCOR.  It  is  structured  in  modular  form  and  is 
based  on  a  keyword  input  format.  This  arrangement  facilitates  making  a  large 
number  of  runs  using  different  data  input,  different  output  formats  and  a 
variety  of  computations  while  requiring  simple  keyword  card  changes.  i 

The  program  reads  data  from  a  permanent  file,  sets  up  arrays  of 
specified  size  and  computes  the  auto  correlation  or  cross  correlation  matrix 
for  the  two  arrays  via  a  very  efficient  two  dimensional  Fast  Fourier  Transform 
(FFT)  code.  The  algorithm  was  between  5  and  20  times  faster  than  other  FFT 
implementations  on  the  benchmark  cases  we  ran.  Output  from  the  use  of  a  two 
dimensional  square  wave  test  case  was  fully  consistent  with  theoretical 
analysis.  This  test  was  used  because  of  its  simple  structure,  and  well 
known  properties  while  at  the  same  time  providing  a  good  check  on  the 
accuracy  and  correctness  of  the  code. 

During  the  development  of  the  WETCOR  package  parallel  efforts  con¬ 
sidered  alternatives  to  the  FFT  approach.  Much  attention  was  given  to  methods 
which  offered  savings  in  processing  time.  Contemporary  journals  on  informa¬ 
tion  and  image  processing  contain  an  abundance  of  references  to  Walsh,  Haar 
and  other  orthogonal  families  of  functions  and  their  corresponding  "fast 
transforms."  The  greatest  strength  of  both  the  Haar  and  Walsh  transforms 
is  the  reduction  in  processing  times  by  factors  of  2  to  20  over  comparable 
times  using  the  FFT  algorithm.  However,  as  related  to  this  specific  problem 
they  are  not  very  useful,  at  the  present  time.  The  principal  drawback  is 
the  difficulty  in  interpreting  and/or  decoding  the  convolution/deconvolution 
equations  which  are  used  in  the  discrete  fast  Haar  and  fast  Walsh  transforms. 

WETCOR  does  correlation  analyses  of  data  arrays  which  are  based  on 
raster  scans  of  weather  patterns.  The  input  data  are  stored  in  bytes,  each 
containing  the  pixels  from  a  raster.  A  portion  of  an  input  raster  is  select¬ 
ed  as  a  base  array.  Another  array  is  selected  and  correlated  against  the 
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base  array  using  the  Fast  Fourier  Transform. 

The  program  prints  the  initial  data  and  the  final  correlation  array 
in  a  compact  exponential  form.  Other  options  for  printing  intermediate 
results  and  for  handling  the  data  differently  are  within  the  program's 
capabilities . 

The  program  structure  is  modular.  Additions  and  modifications  are 
easily  made. 

The  main  structure  of  WETCOR  is  modular  with  the  choice  of  the  order 
of  execution  of  the  modules  controlled  by  the  order  of  the  input  cards. 

Most  of  the  modules  reset  parameters.  One  of  the  main  modules  reads  in  and 
sets  up  the  base  array.  The  other  reads  in  the  second  array  and  does  the 
correlation  using  the  Fast  Fourier  Transform. 

The  size  of  the  array  extracted  from  the  file  WETDAT  is  variable.  This 
array  is  defined  by  its  time  and  the  position  of  its  corner.  If  its  size 
extends  beyond  the  array  stored  on  WETDAT,  zeros  are  filled  in.  This  array  is 
then  packed  into  an  array  whose  dimensions  must  be  powers  of  2,  which  is  a 
restriction  imposed  by  the  Fast  Fourier  Transform  routines  used.  The  average 
is  subtracted  from  the  array  and  the  edges  are  rounded  to  avoid  aliasing 
effects . 

Both  arrays  are  Fourier  transformed;  the  transforms  combined;  and  the 
reverse  transform  performed.  This  final  array  is  printed.  It  may  also  be 
searched  for  the  maximum  and  a  scaled  output  produced  which  emphasizes  the 
maximum.  A  centered  output  is  available.  The  final  array  may  be  output  on  a 
special  file. 
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INPUT 


A  -  Punched  Card 


The  input  to  WETCOR  is  free  field  keyword  input.  Each  card  consists 
of  a  keyword  followed  by  any  number  of  fields.  The  keyword  and  fields  are 
separated  by  break  characters:  either  "("  or  .  The  scanning  of  each 

card  is  terminated  by  either  ")" »  the  end  of  the  card,  or  followed  by 
20  blanks.  The  remainder  of  the  card  can  be  used  for  comments.  (A  card  with 
C  in  the  first  column  is  treated  as  a  comment  card  which  wil 1  appear  in  the 
output) . 

The  basic  deck  consists  of: 

1st  -  title  card  -  80  columns  used  for  labelling. 

2nd  -  Nth  -  keyword  cards  -  (control  the  type  and  amount  of 

processing  to  be  done) . 

N  +  1st  -  End  of  7/8/9  -  (terminates  the  run). 

The  keywords  are  described  below.  The  fields  are  indicated  as  I  for 
integer,  R  for  real,  and  H  for  hollerith  or  character  strings.  An  integer 
has  only  a  sign  and  an  integer  string.  A  real  has  a  and/or  "E".  A 
field  containing  any  other  character  will  be  interpreted  as  a  character 
string.  If  an  input  field  has  no  characters  or  only  blanks,  it  is  treated 
as  not  being  given  and  default  values  are  assigned.  For  numerical  data 
this  usually  means  the  value  zero  is  used.  The  default  values  given  below 
are  the  values  used  if  that  field  is  not  used  or  the  values  used  if  that 
keyword  is  not  used  to  set  the  values.  The  integers  in  parentheses  indicate 
the  field  used;  e.g.  for  the  first  case  given:  BASE  (  ,7311,2261)  does  not 
use  the  first  field;  column  and  line  numbers  are  given  in  the  next  2  fields; 
the  remaining  fields  are  not  used. 
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4m 


The  keywords  available  can  be  divided  into  four  basic  groups: 

1)  Keywords  which  must  be  used  in  all  runs. 

2)  Keywords  which  control/identify  the  structure  of  the  data  being 
used. 

3)  Keywords  which  identify  the  operations  to  be  performed. 

A)  Keywords  which  control  the  output  (format  and  amount  of 

information) . 


1.  Keyword  cards  required  for  every  run:  BASE,  GO 

(a)  BASE  (Al,  A2,  A3,  A4,  A5,  A6,  A7,  A8)  -  the  Base  Keyword  reads 

in  the  first  base  array  and  adjusts  it  to  specifications. 

Al  -  integer  -  time  in  form  hhmm.  If  Al  is  not  specified  the 
next  data  file  is  used. 

If  Al  is  less  than  -2400,  data  is  generated  for 
testing  (see  below) . 

A2  -  integer  -  column  number  ) 

A3  -  integer  -  line  number  (  of  the  upper  left  hand  comer 

of  the  picture 

A4  -  integer  -  data  value  to  be  used  (default:  1000) 

A5  -  integer)-  lower  limit  -  x  axis  (to  be  set) 

A6  -  integer)-  upper  limit  -  x  axis  (to  be  set) 

A7  -  integer)-  lower  limit  -  y  axis  (to  be  set) 

A8  -  integer!-  upper  limit  -  y  axis  (to  be  set) 

(b)  GO  (Al,  A2,  A3,  A4 ,  A5)  -  The  GO  keyword  reads  in  the  second 

array,  adjusts  it  and  does  correlation  calcula¬ 
tions  . 

Al  -  integer  -  time  in  form  hhmm.  If  it  is  not  given,  the  next 
data  file  is  used. 

A2  -  integer  -  column  number 

A3  -  integer  -  line  number 

A4  -  integer  -  column  shift  (if  desired) 

A5  -  integer  -  line  shift  (if  desired) 
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The  shifts  are  applied  to  the  data  portions  on  WETDAT  when  new  data 
is  requested.  The  shifts  are  applied  to  the  stored  data  when  autocorrelation 
is  being  calculated. 


2.  Keyword  cards  which  control/ identify  data  structure. 


EVERY  (Al,  A2,  A3,  A4)  -  indicates  column  and  line  sampling 

rate  on  WETDAT. 


Al  -  integer  -  column  rate 


A2  -  integer  -  line  rate 


of  array  to  be  correlated 
(Go  array) 

(default:  use  old  values) 


A3  -  integer  -  column  rate  I  of  base  array 


A4  -  integer  -  line  rate 


REWIND 


INDIM  (Al,  A2,  A3) 

Al  -  integer  - 

A2  -  integer 
A3  -  integer 

DIM  (Al,  A2 ,  A3) 


(default  if  A3,  or  A4  is  blank: 
use  values  for  other  array. 

Default  if  the  card  is  terminated 
early;  use  old  values) 

The  defaults  at  load  time  are  1, 

2,  1,  2. 

rewinds  the  file  WETDAT. 

Searching  is  only  done  in  forward 
direction;  if  an  earlier  file  is 
to  be  used  a  rewind  is  necessary 
before  searching. 

sets  input  array  dimensions.  This 
is  not  restructed  to  powers  of  2 
but  will  usually  be  1/2  of  DIM 
value . 

set  to  1  for  GO  array 
set  to  2  for  BASE  array 

column  dimension 

line  dimension  (default:  column 
dimension).  The  defaults  at 
load  time  are  all  32. 

sets  dimensions  of  working  arrays. 
These  must  be  powers  of  2.  The 
card  fields  are  identical  to  those 
of  INDIM.  The  defaults  at  load 
time  are  all  64. 


uj. 


(e)  PRP  (Al,  A2,  A3,  AA,  A5,  A6,  A7,  A8)  -  sets  the  parameters 

controlling  the  routine  which 
prepares  the  input  arrays  for  the 
Fast  Fourier  transformation. 

Empty  fields  are  treated  as  zeros. 
The  load  time  defaults  for  both 
arrays  are  given  in  parentheses 
at  the  right. 

Al  -  integer  -  choses  array  as  with  INDIM 

A2  -  integer  -  ,  set  to  0:  center  array 

(set  to  1:  use  shift  values  in  5  and  6 
default:  1 

A3  -  integer  -  edge  treatment. If  several  choices  are  wanted, 

sum  the  corresponding  values. 

The  edges  of  the  input  array  are 
scaled  by  cos^.  This  edge  round¬ 
ing  can  be  eliminated  by  defining 
the  edge  size  to  be  0.  The  edge 
treatment  parameter  values  are: 

1:  Use  values  7  to  define  edge  size  in  x-direction 
instead  of  — jg-  dimension. 

2:  Use  value  8  to  define  edge  size  in  y-direction  instead 
of — jg-  dimension. 

A:  Zero  edge  of  output,  not  the  edge  of  the  minimal  array. 
8:  Do  not  propagate  the  edges  of  the  input  array  (working 
array  set  to  zero  outside  input  range) . 

16:  Set  input  array  edges  equal  to  input  average, 
default:  16 

AA  -  integer  -  Normalization  parameter 

0:  Subtract  average  and  divide  by  average 
1:  Divide  by  average;  do  not  subtract  average 
2:  Subtract  average;  do  not  divide  by  average 
3:  Do  not  normalize 
default:  2 

A5  -  integer  -  z-shift  Used  if  field  2  is  non  zero 

A6  -  integer  -  y-shift 
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A7  -  integer  -  x  edge  size  Used  if  field  3  requires  edge  sizes. 

A8  -  integer  -  y  edge  size 

default:  +  for  A5  through  A8:  0. 


Keyword  Cards  which  identify  the  operations  to  be  performed. 

(a)  PEAK  } control  the  search  for  a  peak  in  the  correlation 
NO  PEAk( array  (default  at  load  time:  NOPEAK) 


(b)  AUTO  )  control  the  generation  of  the  second  array  from  the 
(  base  array 

NAUTO  j  for  autocorrelation  testing  (default  at  load  time: 
NAUTO). 

Keyword  Cards  controlling  output 


(a)  PRINTH  (VI,  V2,  ....  Vn)  -  controls  which  arrays  are  to  be 
printed.  Each  field  is  an  integer  indicating  a 
particular  array  to  be  printed  during  the  calculation. 
The  default  is  (12,  11,  6).  The  acceptable  values  for 
Vj  and  the  corresponding  output  are  as  follows: 


1:  Working  version  of  GO  array 

2:  Working  version  of  Base  array  (printed  for  each 
calculation) 

3:  Transform  of  array  to  be  correlated 
4:  Transform  of  base  array 

5:  Direct  product  of  transforms 

6:  Inverse  transform  of  product;  (correlation  array) 

7:  Correlation  array  shifted  so  that  origin  is  centered 

11:  GO  array  to  be  correlated  as  read  from  WETDAT 

12:  Base  array  as  read  from  WETDAT 
13:  Working  version  of  base  array  when  generated 
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The  output  arrays  are  in  a  compressed  format.  Each  array  entry  is 
multiplied  by  the  scaling  factor  printed  at  the  top  of  the  array.  The 
three  most  significant  difits  greater  than  1  are  then  output  with  a  leading 
symbol  indicating  the  sign  and  the  power  of  10  which  has  been  supressed. 

Symbols  for  positive  blank  +  *  P  P  *** 

negative  -  =  X  M  M  *** 

10  **  0  12  3  4 

(b)  PRINTB  (Al)  Prints  prepared  base  array 

AL  REAL  Scale  factor  to  be  used.  (Default  is  1.) 

(c)  SCALE  (Al,  A2)Sets  output  scale  factors.  These  are  adjusted 

during  calculation. 

Al  -  REAL  -  Scale  for  GO  array  to  be  correlated 

A2  -  REAL  -  Scale  for  BASE  array.  (Default:  same  scale  for 

both) . 

The  load  time  default  for  both  is  1. 


(d) 

TITLE 

Causes  the  next  card  to  be 

read  and  used  as  a 

new  title  card. 

(e) 

PRINT 

I  Causes  the  keyword  cards  to 

printed. 

1 

be  printed  or  not 

NOPRINT 

‘  Default  is  PRINT 

WETDAT  is  read  by  using  unformatted  READ  statements.  The  file  is 
structured  into  records  containing  integer  values  as  follows: 

Record 

1  hour  and  initial  column  number 

2-92  line  number,  data  values  for  that  line 
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Ninety-one  lines  are  expected  -  each  of  which  is  expected  to  have  at 
least  90  values .  These  numbers  are  in  common  /WETDAT/ .  The  column  and  line 
numbers  are  expected  to  change  in  steps  of  1  and  2  respectively,  unless 
changed  by  the  EVERY  keyword.  This  could  be  made  dependent  on  the  more 
recent  WETDAT  files  which  have  the  steps  and  other  information  in  a  longer 
first  record. 

CORARY  is  written  by  unformatted  WRITE  statements.  Each  array  is  a 
file  containing  3  records.  The  first  record  consists  of  120  characters 
containing  the  data  and  time  when  the  file  is  written.  The  second  record 
contains  16  quantities:  the  array  dimensions,  the  sampling  rates,  the 
operative  keyword  (probably  GO) ,  and  an  array  of  11  items  containing  the 
input  fields.  The  third  record  is  the  data  array. 

The  normal  output  from  WETCOR  consists  of  the  keyword  cards  which 
control  the  run  and  the  arrays  specified  by  those  keywords.  The  arrays 
are  given  in  the  abbreviated  numerical  form  described  under  PRINT  H.  The 
correlation  array  may  be  output  on  CORARY. 

A  number  of  methods  exist  which  generalize  the  Fourier  transform 
either  by  considering  the  class  of  Fast  Unitary  Transforms  (which  include 
the  FFT)  or  by  considering  group  characters.  The  practical  interest 
here  is  in  the  computational  efficiency  inherent  in  these  more  general 
isforms. 

For  an  arbitrary  sequence  of  functions,  the  Gram-Schmidt  process 
generates  a  sequence  of  orthogonal  functions.  Any  continuous  function  can 
be  expressed  via  this  orthogonal  set  with  minimum  L ^  error;  for  certain 
classes  of  functions  the  convergence  is  uniform.  Most  important  among 
these  are  the  Haar  functions,  which  assume  2  values ,  and  the  Walsh  functions 
with  values  +  1.  Both  have  discrete  analogies  and  Fast  Discrete  Haar/Walsh 
transforms  are  computable,  FHT  and  FWT.  Because  of  the  simplicity  of  the 
basis  functions  much  greater  computational  savings  can  be  obtained  than  from 
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FFT  (up  to  30  times  faster  than  FFT) . 

Advantageous  use  can  be  made  of  FHT/FWT  in  certain  applications:  date 
transmission/reconstruction  in  which  one  represents  a  given  signal  in  some 
sense  and  reconstruct  it  from  the  minimal  representation.  A  number  of  serious 
difficulties  arise  with  these  transforms  due  to  the  fact  that  the  relation 
with  the  circle  has  been  lost.  For  example:  1)  no  natural  interpretation 
in  terms  of  frequency  exists  (the"sequencyM  viewpoint  of  Harmuth  for  Walsh 
functions  lacks  physical  meaning);  2)  due  to  absence  of  the  circle  relation¬ 
ship  the  important  convolution  theorem  is  not  available  (forced,  analogies 

to  a  convolution  theorem  via  dyadic  convolutions  have  been  made  but  their 

»■ 

interpretations  are  not  clear) . 

The  striking  advantages  of  FWT  and  FHT  over  the  usual  FFT  in 
computational  effort  should  motivate  further  investigation  in  this  area. 

The  historical  situation  regarding  orthogonal  functions  at  the 
beginning  of  the  twentieth  century  was  one  of  well  known  and  useful  kinds 
of  such  function:  the  trigonometric  functions  which  occur  in  Fourier  series; 
orthogonal  polynomials  such  as  those  of  Legendre,  Hermite,  and  Laguerre; 
Bessel's  functions,  the  Sturm-Liouville  series,  and  other  special  functions. 
But,  there  was  no  general  theory  embracing  all  such  systems  of  functions. 
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